Spatial Direct Numerical Simulation of Boundary-Layer 
Transition Mechanisms: Validation of PSE Theory 

R.D. Joslin and C.L. Streett 

NASA Langley Research Center, Hampton, VA 23681 

and 

C.-L. Chang 

High Technology Corporation, Hampton, VA 23666 

Abstract. A study of instabilities in incompressible boundary-layer flow on a flat plate 
is conducted by spatial direct numerical simulation (DNS) of the Navier-Stokes equations. 
Here, the DNS results are used to critically evaluate the results obtained using parabo- 
lized stability equations (PSE) theory and to study mechanisms associated with breakdown 
from laminar to turbulent flow. Three test cases are considered: two-dimensional Tollmien- 
Schlichting wave propagation, subharmonic instability breakdown, and oblique-wave break- 
down. The instability modes predicted by PSE theory are in good quantitative agreement 
with the DNS results, except a small discrepancy is evident in the mean-flow distortion 
component of the 2-D test problem. This discrepancy is attributed to far-field boundary- 
condition differences. Both DNS and PSE theory results show several modal discrepancies 
when compared with the experiments of subharmonic breakdown. Computations that allow 
for a small adverse pressure gradient in the basic flow and a variation of the disturbance 
frequency result in better agreement with the experiments. 
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1. Introduction 


In the past century, numerous attempts have been made to understand and predict the 
transition from laminar to turbulent flow in boundary layers. Most of this effort stems from 
the independent theoretical accomplishments of Orr (1907) and Sommerfeld (1908) at the 
turn of the 20th century. The Orr- Sommerfeld equation is based on linearized disturbance 
equations and is a successful example of classical hydrodynamic stability theory. Nearly 20 
years later Tollmien (1929) solved the Orr- Sommerfeld equation, which led to the calcula- 
tion of a critical Reynolds number for the onset of instability. On the same road, Schlichting 
(1933) computed amplification rates of disturbances in the boundary layer. The first exper- 
imental confirmation of the theory was accomplished by Schubauer and Skramstad (1943, 
1947); a vibrating ribbon was used to impress a disturbance into the boundary layer, and 
hot wires were used for measurements. These contributions (and others) span over some 40 
years of research, after which theory and experiments finally agree on the initial growth of 
disturbances in boundary layers. 

Stability theory has gained wide acceptance and is now a well-established tool in the 
research and engineering community. The first reasonably comprehensive method for pre- 
dicting transition was derived from stability theory: the e N method by Smith and Gamberoni 
(1956) and van Ingen (1956). Although the e N method is widely used to predict transition 
in a broad class of flows, this method has considerable shortcomings. Because the approach 
uses linear stability theory, the limitations of the theory are imbedded in the methodology. 
Specifically, a quasi-parallel boundary layer is assumed, and because the equations are linear, 
no amplitude information about the ingested disturbance can be taken into account. Finally, 
the method is semi-empirical; therefore, some foreknowledge is required of the flow in tran- 
sition. The true physical problem involves disturbances that interact in a nonlinear manner 
in later stages of transition and are imbedded in a growing boundary layer. Consequently, a 
method that accounts for nonparallel flow and nonlinear interactions is necessary to predict 
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transition. At present, such an all-encompassing method of transition prediction is beyond 
our grasp; however, progress has been made in recent years. 

In the last decade, significant accomplishments have been made in the theoretical pre- 
diction of transition beyond the linear growth stage. Based, in part, on pioneering attempts 
at nonlinear theories (e.g., Benney and Lin, 1960; Craik, 1971), Orszag and Patera (1980, 
1983) by numerical experiments and Herbert (1984) by theory discovered a secondary insta- 
bility mechanism. Herbert’s (1984) theory is based on Floquet theory and accounts for an 
experimentally observed three-dimensional (3-D) instability. Although the governing equa- 
tions are linearized and a local parallel-flow assumption is made, remarkable agreement is 
obtained between predictions from Herbert’s theory and experimental results, particularly 
for the peak-valley splitting mode identified by Klebanoff et al. (1959, 1962) and for the 
peak-valley alignment mode observed by Kachanov et al. (1977, 1984). These modes are 
two distinct and different routes to transition that are discriminated based on the initial 
disturbance levels. After nearly a decade of verification, Herbert’s theory for secondary in- 
stabilities is generally accepted and widely used by the research community to understand 
certain aspects of transition in boundary layers. 

More recently, Herbert and Bertolotti (1987) have devised a nonlinear, nonparallel com- 
putational method based on so-called parabolized stability equations (PSE). The benefits 
and limitations of this new theory have yet to be determined and will be explored somewhat 
in the present paper. Before the development of this theory, the only approach to solve 
the nonparallel, nonlinear boundary-layer transition problem was by direct numerical sim- 
ulation (DNS), except in the large Reynolds number limit where asymptotic methods have 
had some success (e.g., Smith, 1979; Hall and Smith, 1991). To date, most studies using 
DNS have been limited to the temporal formulation-where a spatially-periodic computa- 
tional domain travels with the disturbance and the temporal evolution of the disturbance is 
computed. This enabled simulations into the later stages of transition (Zang and Hussaini, 
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1987,1990; and Laurien and Kleiser, 1989), thus providing a database of qualitative infor- 
mation, which, however, lacks the physically-realistic spatial representation. Spatial DNS 
computes spatially-evolving disturbances and would provide needed quantitative information 
about transition. Progress in spatial DNS has been made by, among others, Danabasoglu et 
al. (1990, 1991) for channel flow and Fasel (1976), Spalart (1989), Fasel et al. (1990), Rai 
and Moin (1991a, 1991b), and Joslin et al. (1992a) for boundary-layer flow. Results obtained 
from spatial DNS have been compared qualitatively and, with some success, quantitatively 
to linear stability theory, secondary instability theory, and available experimental results. 
For a more complete list of accomplishments in transition prediction with DNS, refer to the 
review by Kleiser and Zang (1991). 

The goal of the present research is to quantify various transitional-flow mechanisms of 
interest from the linear region to the nonlinear breakdown stage and to provide a critical 
comparison of spatial DNS and PSE theory results. To date, this comparison offers the most 
rigorous test of the PSE approach for accuracy; the main focus will be to point out strengths 
and potential weaknesses of PSE theory and the impact of these weaknesses on the overall 
flow-field prediction. To accomplish this goal, three test cases are computed by spatial DNS 
and then compared to PSE theory: two-dimensional (2-D) Tollmien-Schlichting (TS) wave 
propagation, subharmonic breakdown, and oblique-wave breakdown. 

2. Governing Equations 

To compute the disturbance development, the incompressible Navier- Stokes equations 
are solved. The streamwise direction is x, the wall-normal direction is y, and the spanwise 
direction is z. Instantaneous velocities u = (u, 7, w) and pressure p are decomposed into the 
basic components U_ = (£7, V, W) and P and the disturbance components u = (it, v } w) and 
p so that 

u(x, t) = U_{x_) + u(x, t) and p(x_, t) = P(x )+p(x,t) (1) 

where x_ = (x, y, z) and t is time. 
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The basic flow on a flat plate is generally governed by the full Navier- Stokes equations. 
However, according to order-of-magnitude analysis, the streamwise parabolized form of the 
equations (boundary-layer equations) is adequate. A marching algorithm may be used to 
solve the equations for the basic flow (H, V), or the widely used Blasius similarity profile can 
be employed. These approaches have been verified by comparisons with experiments (see 
Schlichting, 1955). For the present study, the Blasius flow is used for the basic state. 

To determine the disturbance component of equation (1), the equations that are solved 
result from substituting equation (1) into the Navier-Stokes equations and subtracting the 
basic flow equations. These unsteady, nonlinear disturbance equations are 
Q u 1 

-= + (u • V)u + (U • V)u + (u • V)U = -Vp + — ' V 2 u (2) 

ut 

and the continuity equation 

V-u = 0 (3) 

The equations are nondimensionalized with respect to the freestream velocity U 0 Q , the kine- 
matic viscosity i/, and some length scale at the inflow (e.g., displacement thickness 6*). A 
Reynolds number can then be defined R$* = Uqo8*/u. 

At the wall and in the far field, boundary conditions for equations (2) and (3) are 

u = 0 at y = 0 and u — > 0 as y — > oo (4) 

To prevent spurious reflections, outflow conditions are provided by the buffer-domain tech- 
nique (Streett and Macaraeg, 1989; Joslin et al. 1992a). For the present study, the dis- 
turbance forcing takes the form of eigenfunctions that are imposed at the inflow boundary. 
Controlled disturbances at the inflow are required to facilitate the comparison. To this end, 
the inflow condition for both DNS and PSE is given by the basic flow plus the disturbance 
forcing functions 

N z N t 

U 0 = U 0 + A °m,n-[M-m,n(y) eX (8 z ~ mut)]] (5) 

n= — N z m= — N t 
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where U_ 0 is the inflow basic component; A ° m n are prescribed 2-D and 3-D disturbance am- 
plitudes; f3 is an imposed spanwise wave number; and u is an imposed disturbance frequency. 
The terms n (y) are complex eigenfunctions found either by solving the Orr-Sommerfeld 
and Squire equations or obtained from a secondary instability theory (see Herbert, 1984). 
The eigenfunctions u ° m n (y) are normalized with respect to the maximum streamwise- velocity 
component such that the initial amplitudes of the induced disturbances are prescribed by 
4° 

^ 171,71 * 

3. Computational Methods of Solution 

In this paper, spatial DNS is used to study breakdown mechanisms occurring in tran- 
sitional flows, and to evaluate the accuracy of PSE theory in predicting the evolution of 
convective disturbances. A description of both the DNS and PSE theory approaches follows. 


Spatial DNS. 

In a previous study by the authors described in Joslin et al. (1992a), the spatial DNS 
approach involved fourth-order compact differences used on a nonstaggered grid. For this 
study, sixth-order compact-difference schemes have replaced the fourth-order schemes, and 
a staggered grid has replaced the nonstaggered-grid formulation. 

In the streamwise direction (^-direction), fourth-order central finite differences are used 
for the pressure equation. At boundary and near-boundary nodes, fourth-order differences 
are used. For first and second derivatives in the momentum equations, sixth-order compact 
differences are used. As described by Lele (1992), the difference equations are 

-/)_i + fl + -f'i+i = ~ -^-(fi+i ~ fi-i) + gg^ (/*+ 2 - fi — 2 ) (6) 


and 


jjfi-i + /; 


+ iP‘+i = 


12 

Tuf 


(/*+ 1 — 2 f z + fi-l) + 2 (fi+2 — 2 ft + fi — 2 ) (7) 


44 hi 


where h x is the uniform streamwise step size, and / is an arbitrary function whose derivatives 
are sought. At boundary and near-boundary nodes, explicit fifth-order finite differences are 


6 



used. Equations (6) and (7) lead to tridiagonal systems which can be solved efficiently by 
LU-decomposition with appropriate backward and forward substitutions. 

In the wall-normal direction (y-direction), Chebyshev series are used to approximate the 
disturbance at Gauss-Lobatto collocation points. Since this series and its associated spectral 
operators are defined on [-1, 1] and the physical problem of interest has a semi-infinite domain 
[0, oo) or a truncated domain [0,y max ], a transformation is employed. Studies of spectral 
methods and mapping tranformations in unbounded regions have been conducted by Grosch 
and Orszag (1977) and Boyd(1989). Here an algebraic mapping is used: 


VmaxSp i 1 + y) 

2 s p + y max (! - y) 


( 2 Sp T Umax^y y 

or y = 7 r 

t/max \Sp "by) 


maxSp 


(8) 


where y £ [0, y max ) and y £ [—1, 1]; y ma x is the wall-normal distance from the wall to the far- 
field boundary in the truncated domain, and s p controls the grid stretching in the wall-normal 
direction. In the spanwise direction (^-direction), periodicity is assumed, allowing for Fourier 
series representations. Using Fourier series, spectral accuracy is obtained in the spanwise 
direction and Fast Fourier Transforms (FFT) or sine/cosine transforms may be used, allowing 
for the fast computing of derivatives. For more details on the spectral methods used here refer 
to Canuto et al. (1988) and Joslin et al. (1992a). Although spectral accuracy is obtained in 
the wall-normal and spanwise directions, the overall spatial-discretization scheme is formally 
fourth-order accurate as a result of the streamwise discretization. 

For time marching, a time-splitting procedure was used with implicit second-order 
Crank-Nicolson differencing for normal diffusion terms; an explicit third-order three-stage 
Runge-Kutta (RK) method was used for the remaining terms. This time-stepping procedure 
was used successfully by Streett and Hussaini (1991) for Taylor- Couette flow simulations. 
The pressure is omitted from the momentum equations (2) for the fractional RK stage, 
leading to 

?/ * _ q .m f<m 

- - = CRH m (u) + ^D 2 (u* + u m ), (9a) 

o 
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where 


H m (u ) = L(u) m + C^H 171 - 1 ^), 


( 96 ) 


and operator 

L(u) = (U_ • V)u + (u • V)U + (u ■ V)u — — — V^ z u, (9c) 

R s * 

Here u* are disturbance velocities at the intermediate RK stages, u m are velocities at previous 
RK stages ( m = 1, 2 or 3), u° are velocities at the previous time-step, = d 2 / dx 2 -\-d 2 / dz 2 , 
h t is the time-step size, and D is the wall-normal spectral derivative operator. For a full RK 
stage, the momentum equations with the pressure are 


u m+1 - u m 


= C?H m (u) + -2- D 2 (u m+1 +u m ) - Vp m+1 
Rs* 


( 10 ) 


Subtracting (9) from (10) leaves 


u m+1 - u * 

K 


= — Vp 


m+1 


(ii) 


By taking the divergence of (11) and imposing zero-divergence of the flow field at each RK 
stage, a pressure equation is obtained 


vV " +1 = ^ r(V '- ) (12) 

which is subject to homogeneous Neumann boundary conditions. This boundary condition is 
justified in the context of a time-splitting scheme as discussed by Streett and Hussaini (1991). 
The solution procedure is as follows: The intermediate RK velocities (u*) are determined by 
solving equation (9). The pressure correction (p m+1 ) is found by solving (12). Then, the full 
RK stage velocities ( u m+1 ) are obtained from (11). Upon solving the above system three 
consecutive times, full time-step (n + 1) velocities are determined, where u n+1 = u m=3 . The 
RK coefficients and time-steps given by Williamson (1980) are 

(Cl Cl Cl\ ( 11/2 0 \ f h]] ( 1/3 h t ] 

Cl Cl Cl = 9/4 1/2 -4 and l h] } = 5/12 h t > (13) 

C 2 3 Cl) V 32 / 15 1/2 -153/32/ (hi J ( l/4h t J 



where the sum of the three RK time-steps equals the full time-step (h t ). 

A significant improvement in the spatial DNS approach over that used in Joslin et al. 
(1992) involves the use of staggered grids in the wall-normal direction. The velocities and 
the 2-D pressure component are computed at Gauss-Lobatto points as before, and the 3- 
D pressure components are computed at Gauss points. With this staggered-grid pressure 
solver, reductions of 10-50 percent in CPU time, 30-50 percent in virtual memory, and 50-70 
percent in run-time disk requirements were obtained, compared with the nonstaggered-grid 
solver. This staggered grid approach affects only the pressure equation (12) and is described 
in the remaining portion of this subsection. 

To obtain the pressure correction (p) for the 2-D and 3-D boundary-layer problems, 
solutions of Poisson equations for each RK stage are required. For 3-D simulations with 
spanwise periodicity assumed, the pressure correction is determined in transform space where 
the Fourier coefficients are evaluated. In transform space, the Poisson equations become 
Helmholtz equations. In order to solve the equations efficiently, a fast elliptic solver is 
required. For this purpose, the tensor-product method described by Lynch et al. (1964) is 
used. On a nonstaggered grid, this approach was employed by Danabasoglu et al. (1990, 
1991) for the channel problem and by Joslin et al. (1992a) for the boundary-layer problem. 

For the present staggered-grid formulation, the discretized equations for the 3-D com- 
ponents are 

(H - l3 2 J) Pn + Pn X T = R n (14) 

where matrix operations are ordered with respect to (y, x). p n are Fourier coefficients of the 
desired pressure solution p, f3 n = n/3 are spanwise wave-number coefficients of the Fourier 
series, I is the identity matrix, X T is the transpose of the streamwise central finite-difference 
operator, and R n are known Fourier coefficients of the right side of the pressure equation 
(12). The following matrix operations determine the wall-normal operator H: 

H = I^DDI^ ( 15 ) 
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where D is a spectral wall-normal derivative operator for the stretched grid, and D is the 
derivative matrix with the first and last rows set to zero. This modification enforces the 
homogeneous Neumann boundary conditions required for p on the Gauss-Lobatto grid. The 
interpolation matrix I° L operates on variables at Gauss-Lobatto points and transforms them 
to Gauss points; the interpolation matrix I° L performs the inverse operation. These opera- 
tors are defined as 


Iol = (C^kkiCo,)^ and = (C~k)j,k(C a )k 


(16a) 


where 


(C G )k,j = = cos [ 


2 r {j + \kk v KJ , 2 , 

(°g )j,k = cos [- 


-O' + \k k ^ 


N yCk 


N , 


N , 


(166) 


with j = 0, 1, . . . , N y — 1; A: = 0, 1, , N y ; and 


{C GL )k,j = = cos [=H, {C GL ) h k = cos 


NyCkCj N t 


N, 


(16c) 


with j = 0,1,..., Ny] and k = 0, 1, ... , N y . The constants ci are defined as ci = 1; i = 
1,2,..., — 1; and c 0 = cj^ y = 2. 

Using the tensor-product method described by Lynch et al. (1964), the matrix H may 
be decomposed into 

H = QAQ- 1 (17) 

where A is a diagonal matrix of eigenvalues and Q is the corresponding matrix of eigenvectors. 
Temporary matrices are introduced and defined: 


Pn — Q Pn and G n — Q R 


-l 


(18) 


Substituting (17) and (18) into (14), one obtains 


(A - 0i)p n + p n X T = G n , 


(19) 


Equation (19) is used to solve for p n , which is then used in (18) to solve for p n . The 
solution is then inverse tranformed to physical space. The operators: 77, Q,Q _1 ,A and 
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X T are all mesh dependent matrices and need be calculated but once, as is the influence 
matrix. Since X T is a fourth-order accurate pentadiagonal matrix, the LU-decomposition 
method is used and provides and efficient method to solve for the pressure. The influence- 
matrix technique is also used for the pressure solver to ensure that the continuity equation 
is discretely satisfied. Details of the influence-matrix technique are given by Streett and 
Hussaini (1991), Danabasoglu et al. (1991), and Joslin et al. (1992a). Since the pressure 
equation (12) is an inviscid calculation, it is well posed provided that boundary conditions on 
the normal component of velocity only are enforced. At the end of each full RK time-step, a 


non-zero tangential velocity component may therefore arise at the computational boundary. 
This is referred to as a “slip-velocity.” To alleviate this problem, intermediate boundary 
conditions as described by Streett and Hussaini (1991) and Joslin et al. (1992a) are used 


and given by 


(1 + 7T^r) V ^ - 


where u BC = 0 for a rigid wall and y LBC = u 0 for an inflow condition or a wall slot condition 


evaluated at the appropriate time in the RK stage. 


Disturbances are introduced into the boundary layer by forcing at the inflow boundary. 
At the outflow, the buffer-domain technique of Streett and Macaraeg (1989) is used. Briefly, 
the streamwise convective terms are attenuated to zero in a small region that is appended to 
the end of the physical domain. As shown by Joslin et al. (1992a) for flat-plate boundary- 
layer flow, a buffer-domain length of approximately three streamwise wavelengths is adequate 
to provide an attenuation function that is smooth enough to avoid the reflection of waves. 
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PSE Theory. 


In contrast to the spatial DNS approach described above in which the complete Navier- 
Stokes equations are solved, the PSE approach seeks approximate solutions to the unsteady 
Navier-Stokes equations, which can be obtained by an efficient marching procedure. Gen- 
erally, there are a number of ways to parabolize the Navier-Stokes equations; however, any 
acceptable approximation must be able to capture the physics of instability waves. One such 
method was first suggested by Herbert and Bertolotti (1987); more details of this method and 
a comprehensive discussion concerning nonparallel and nonlinear effects on Blasius boundary 
layers was given by Bertolotti et al. (1992). The underlying notion of the PSE approach is 
to first decompose the disturbance into an oscillatory wave part and a shape function part. 
By properly choosing a streamwise wave number to resolve the wave motion, the govern- 
ing equations reduce to a set of partial differential equations for the shape functions, which 
vary slowly in the steamwise direction and their second-derivatives are assumed negligible. 
These PDE’s are then parabolized by neglecting the dependence of convected disturbances on 
downstream events and by neglecting the second derivatives (d 2 / dx 2 ) of the shape functions. 
Since most of the oscillatory wave motion is absorbed in the streamwise wave number and 
the terms neglected in the shape function equations are of order 1 / R 2 , the resulting system 
should yield the desired physical results. A brief discussion of the theory is as follows. 

For disturbances that are present in the flow field, periodicity is assumed both in time 
and in the spanwise direction. The total disturbance can then be described by the following 
Fourier-series expansion: 

N z N t 

{u,p}(x,y,z,t) = E E {u m!n ,p m ,n}(x,y)exp[t(n/3z - mtot)] (21) 

n= — N z m= — N t 

where N z and N t are the numbers of modes retained in the truncated series, u is an imposed 
frequency, and f3 is an imposed spanwise wave number. Equation (21) is substituted into the 
governing equations (2) and (3), a set of elliptic equations are obtained for the transformed 
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variables {u_ m n ,Pm,n}- Because of the wave nature of these transformed variables, a further 
decomposition is made into a fast-oscillatory wave part and a slow-varying shape function: 

Pm,n\ = {u. mjn ■, Pm,n } exp [z I Oi m ^ n dx\ (22) 

J x 0 

In equation (22), the fast-scale variation along the streamwise direction x is represented 
by the streamwise wave number a mjn ; therefore, the second-order variation of the shape 
function in x is negligible. This observation leads to the parabolized stability equations for 
the shape functions {u m n} p m , n }- These equations are obtained by neglecting all second 
derivatives in the streamwise direction and the terms associated with upstream influence. If 
a njTn are chosen properly, the evolution of disturbances can be described by the parabolized 
equations for the shape functions. Based on decompositions (21) and (22), linear PSE can 
be derived for any disturbance with a given frequency and spanwise wave number. 

To solve the nonlinear problems, the nonlinear convection terms are placed on the right 
side of the governing equations: 


F_(x, y, z, t) = (u ■ V)u. (23) 

For the PSE approach, the governing equations are solved in wave number space. The 
Fourier coefficients, which are obtained from the corresponding Fourier transform of F_ in 
equation (23), provide a nonlinear forcing to each of the linearized shape-function equations. 
These inhomogeneous equations for the shape functions are solved by applying a marching 
procedure along the streamwise direction for each Fourier mode. 

Both Bertolotti (1991) and Chang et al. (1991) have extended the PSE numerical 
approach to the study of compressible boundary layers. The PSE code developed by Chang 
et al. (1991) with M 0 Q ~ 0 is used for the present comparison with both the DNS results and 
the experimental results. A second-order backward differencing is employed to integrate the 
equations in the streamwise direction, and fourth-order, finite-difference schemes are used to 
discretize the normal derivatives. 
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DNS and PSE Boundary Condition Treatment. 


For the DNS approach, the Navier-Stokes equations with imposed boundary conditions 
are solved for full disturbance components, which is equivalent to the left side of equation 
(21). For PSE theory, equations (21) and (22) are substituted into the Navier-Stokes equa- 
tions (the DNS equations) to produce a coupled system of equations for the unknown Fourier 
coefficients (h m , n ,Pm,n)- Boundary conditions are imposed for each Fourier-coefficient equa- 
tion. Because the PSE theory solves a system of equations in Fourier-coefficient space, more 
degrees of freedom exist compared with the DNS. 

At the wall, no-slip boundary conditions (i.e. , homogeneous Dirichlet conditions) are 
enforced for the disturbance equations used by the DNS and the Fourier-coefficient equations 
used by PSE theory. So, both DNS and PSE incorporate the same boundary conditions at 
the wall. In the far field, homogeneous Dirichlet boundary conditions are imposed for the 
DNS computations. This far-held condition is exact at infinity, but to computationally solve 
the system using DNS the semi-infinite domain is truncated. 

For the PSE approach, homogeneous Dirichlet boundary conditions are used for all 
Fourier-coefficient equations except for the mean-how distortion equations. This nonzero 
component arises from assumptions of PSE theory. Unlike the DNS which solves the full 
Navier-Stokes equations, PSE theory reduces the equations to a simplihed parabolic system 
in a manner described in the previous subsection. As a result of this PSE simplihcation, 
the mean-how distortion equations become essentially of the boundary-layer equation type. 
With boundary-layer equations, the wall-normal velocity component approaches a constant 
in the far-held. Similarly, the mean-how distortion equation in PSE theory, which is the 
boundary-layer equation type, incorporates a Neumann boundary condition for the wall- 
normal velocity component. This Neumann condition allows the total normal velocity (mean 
how + mean-how correction) in the far-held predicted by PSE theory to vary at inhnity. 

Thus, the far-held boundary conditions used by both the DNS and PSE approaches 
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are approximate. The spatial DNS far-held boundary conditions cannot be changed to 
mimic the PSE approach because the DNS cannot accomodate an a priori Fourier-modal 
boundary- condition treatment as is present in PSE theory. Changing the Neumann condition 
to Dirichlet resulted in numerical instabilities in the PSE approach. As a result, the far- 
held boundary conditions for PSE theory are different from the boundary conditions for the 
DNS approach for the mean-how distortion equation. With the present boundary-conditions, 
the PSE theory approximate far-held conditions should prove more accurate for a far-held 
boundary hxed close to the wall, and the DNS conditions should prove more realistic for the 
boundary far from the wall. 


4. Results 

The spatial DNS approach has been tested against linear stability theory for grid and 
time-step requirements as described by Joslin et al. (1992a). Excellent agreement in the 
amplitudes and phase of disturbances was found between the DNS results and linear theory; 
those results will not be repeated in this paper. Rather, three vastly different instability 
cases have been selected for study: two-dimensional Tollmien-Schlichting wave propagation, 
subharmonic breakdown, and oblique-wave breakdown. The reasons for the selection of each 
case will be discussed in each subsection below. These test cases serve to provide a com- 
prehensive evaluation of the effectiveness of PSE theory in predicting convective disturbance 
development in boundary-layer flows on flat plates. 

2- 1 ) Tollmien-Schlichting Wave Propagation. 

For this test case, the evolution of an ingested 2-D Tollmien-Schlichting wave is predicted 
by PSE theory and compared with the DNS results. Because the nonlinear evolution of a 
single 2-D wave offers a less complex instability mechanism for the boundary-layer problem, 
this test case is a logical first choice for evaluating PSE theory. 

Bertolotti et al. (1992) performed a similar comparison between DNS and PSE theory 
for the u-velocity components of the fundamental and first harmonic waves. In their work, 
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a body-force term was used in the governing equations to introduce the disturbance for the 
DNS. Then a prescribed inflow (generated by weakly nonlinear theory) was used for distur- 
bance initiation for the PSE. Their results showed good agreement between the fundamental 
wave and first harmonic amplitudes with downstream distance. If a discrepancy in the results 
had occurred, one might attribute the difference in results to either a PSE theory flaw or to 
the different means by- which the disturbances were forced. The intention here is to remove 
the possibility of variations in disturbance forcing, so the results can be compared without 
this potential source of confusion. For the present study, the exact same inflow profiles and 
amplitudes for disturbance forcing are used for both PSE and DNS approaches. As a result, 
the outcome of the comparison should not depend on a variation of disturbance forcing. 

A 2-D Tollmien-Schlichting disturbance with a root-mean-squared (r.m.s.) amplitude 
A\ 0 = 0.0025 is introduced into the boundary layer by a forcing at the inflow for the DNS 
as well as for the PSE calculations. (Note: the definition of n is given in equation (5).) 
Through nonlinear interactions, all other harmonic waves including the mean-flow distortion 
are generated for both DNS and PSE. Calculations are made with an inflow Reynolds number 
R$* = 688.315 and frequency F = 86. To generate resolved benchmark data to test PSE 
theory, the spatial DNS was computed on a grid of 2041 uniformly spaced streamwise nodes 
(60 nodes per disturbance wavelength) and 81 wall-normal collocation points. The outflow 
boundary is 442 8* from the inflow boundary, and the far-field (or free-stream) boundary is 
75 8* from the wall. The DNS parameters were chosen based on convergence studies by Joslin 
et al. (1992a). For the time-marching scheme, the disturbance period is divided into 320 
time steps. For the PSE computational approach, several numerical experiments have been 
performed by varying the grid, far-field boundary location, and the number of Fourier modes. 
These numerical experiments led to the choice of 100 wall-normal grid points; 5 frequency 
modes of series (21) (N t = 6); and a far-field boundary located 58<5* from the wall. 

Figure 1 shows the maximum streamwise amplitudes for the mean-flow distortion rt 0 , 
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fundamental wave u i, and first harmonic 112 predicted by PSE theory and compared to 
the DNS results with the downstream distance. Both the fundamental waves and the first 
harmonics are in good quantitative agreement throughout the initial linear region and the 
later weakly nonlinear region. The PSE results capture early evidence of the first harmonic 
between R$* = 690 and R$* = 900 (barely visible), which is in good agreement with the DNS 
results. The mean-flow distortion components are in good agreement throughout the initial 
linear region; however, a discrepancy begins to occur downstream at an apparent “notch” in 
the results at R$* = 1400. This notch will be explained by looking at the profiles of the mean- 
flow distortion. At the local streamwise location Rs* = 1519, figure 2 shows comparisons 
of the streamwise velocity component. The fundamental (TS) wave and harmonics (not 
shown) are in good quantitative agreement, even in regions of high gradients. Similar to the 
discrepancy in the mean-flow distortion amplitudes shown in figure 1, a comparison of the 
streamwise velocity (u 0 ) profiles illustrates a comparable difference. Although not shown 
here, the comparison of the mean-flow distortion at smaller Reynolds numbers gives much 
better agreement, as expected from figure 1. The notch in the mean-flow distortion results 
identified in figure 1 arises from the change of profile contributions. Initially, the positive 
value of the profile near the wall has a larger amplitude than the retarded, or negative, 
profile amplitude. At the downstream location near R$* = 1400, the negative portion of the 
mean-flow distortion becomes the dominant amplitude, and the amplitude of the positive 
profile near the wall begins to decay. Shown in figure 3, comparisons of the wall-normal 
velocity ( v 0 ) profiles indicate a larger discrepancy between DNS and PSE theory results. 
This discrepancy is due to the homogenous Neumann boundary conditions used in the far 
field for the mean-flow distortion equations in PSE theory. As in the traditional boundary- 
layer equations approach, this boundary condition leads to a nonzero, wall-normal mean-flow 
velocity component in the far-field, as discussed in some detail in the previous section. As 
shown in figure 3, the wall-normal component v 0 approaches a constant for PSE theory. In 
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the present case, the modified far-field normal velocity (i.e., the Blasius solution plus the 
mean-flow distortion) in PSE theory is shown in figure 4 with the Blasius solution. The 
maximum difference in the mean flows occurs near the location of the maximum amplitude 
of the fundamental wave and corresponds to the largest difference in the mean-flow distortion 
components shown in figure 1. 

In this first test problem of 2-D Tollmien-Schlichting wave propagation, the results from 
DNS and PSE theory agree very well for the fundamental and harmonic waves. However, 
a discrepancy exists in the mean-flow distortion component. This discrepancy is a result of 
the difference in the far-field boundary-condition treatment for the two appoaches. 

Subharmonic Breakdown. 

A well-understood breakdown scenario in an incompressible boundary layer on a flat 
plate involves 2-D disturbances and the emergence of aligned and staggered 3-D distinct vor- 
tex structures through spanwise vortex stretching and tilting, or some indistinct, nonunique 
combination of vorticities in the later stages. These vortex patterns are referred to as funda- 
mental, subharmonic, and combination resonant modes, respectively, and may be described 
by secondary instability theory. The goal here is to use this test case for PSE validation; 
however, secondary instability theory will not be discussed. For a good review of secondary 
instabilities, refer to the work by Herbert (1988). For the present study, the subharmonic 
mode of secondary instability breakdown will be computed and the results compared with 
the experiments of Kachanov and Levchenko (1984) (hereafter, K&L). Previously, Fasel et 
al. (1990) used spatial DNS and Herbert (1991) used PSE theory and compared their com- 
puted subharmonic breakdown results with the K&L experimental data. The results for 
both the computation and theory were in relatively good agreement with the experiments. 
In Herbert’s work, comparisons between PSE theory and the spatial DNS results of Fasel et 
al. (1990) agreed well quantitatively for the modes shown; however, as discussed with the 
previous 2-D test case, the PSE and DNS approaches used different methods to introduce 
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the disturbance. The 3-D disturbances were introduced into the DNS computations by wall 
suction and blowing while the PSE disturbances were introduced at the inflow. 

In the present study, the subharmonic breakdown mechanism is revisited. So that a 
detailed evaluation of the PSE theory results is possible, both DNS and PSE theory must 
use the same initial disturbance amplitudes and profiles, which are obtained from linear and 
secondary instability theory. Furthermore, physical features of the K&L experiments differed 
from the earlier computations by both Fasel and Herbert with respect to some of the modal 
trends. In the present paper, an explanation for these discrepancies will be sought. 

For spatial DNS, computations are performed on a grid of 1021 uniformly spaced stream- 
wise nodes, 81 wall-normal collocation points, and 5 symmetric-spanwise nodes. In the 
streamwise direction, the outflow boundary is 442 8* from the inflow boundary; the far field 
boundary is 75 8* from the wall; and the spanwise boundary consists of a length equal to 
one half of the spanwise wavelength, or \ z /2 = 7 r//3. (Note that the spanwise computational 
length would be for the general, nonsymmetric computation.) For the time-marching 
scheme, the disturbance period is divided into 320 time steps, and time is advanced using 
a three-stage Runge-Kutta method. For the PSE computational approach, 100 wall-normal 
grid points are used; 7 frequency modes and 3 spanwise modes from series (21) are used, and 
the far-field boundary is 58<5* from the wall. 

The prescribed primary and subharmonic disturbances are obtained at the Reynolds 
number R$* = 732.711 and the primary frequency F = 124 (which correspond to the exper- 
iments of K&L). The primary wave has an inflow amplitude of A\ 0 = 0.0048 (r.m.s.). The 
subharmonic mode has an inflow amplitude of A° 1 = 0.145 X 10 -4 (r.m.s.) and corresponds 
to a mode with spanwise wave number f3 = 0.2418. 

Figure 5 compares the maximum amplitudes predicted by PSE theory with the results 
from DNS. For this test case, note the extremely good quantitative agreement between PSE 
theory and DNS for the growth rates of the fundamental and the dominant harmonic modes; 
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even the mean-flow distortion components are in good agreement. PSE theory and DNS 
profiles are compared in figure 6 at the local (downstream) Reynolds number R$* = 1067. 
The predictions made by PSE theory are in good overall quantitative agreement with DNS 
results. From these comparisons, the accuracy of the PSE predictions has been established 
for this test case. 

The second goal is to understand the modal discrepancies identified from the comparison 
with the experiments. In figure 7, the DNS results are compared with results of the K&L 
experiments. Furthermore, the computed profiles are in good agreement with the K&L 
experiments as shown in figure 8 at the downstream location where R$* = 1049. Both the 
earlier DNS results of Fasel et al. (1990) and PSE results of Herbert (1991) indicated good 
agreement with the K&L experiments, similar to the present results. But, if the results 
are examined more carefully on a mode-by-mode basis, modal inconsistencies can be found 
and lead to distinct differences in the interpretation of the computational and experimental 
results. In particular, the primary wave (F) and its harmonics (2 F, 3 F) begin to decay much 
sooner downstream in the experiments (Rs* ~ 990) than in the computations (Rs* ~ 1070). 
This observation can be seen clearly in figure 9; the fundamental waves of DNS and the 
experiments are shown with downstream development. A reason for this difference will be 
briefly explored. 

Based on linear stability theory, the discrepancy in the 2-D fundamental mode (shown 
in figure 9) can be attributed to the fact that each wave may have a different frequency. The 
fundamental wave in the experiments reaches a peak amplitude and then decays upstream 
of the theoretical wave, which suggests that the experiment has a larger effective frequency 
than in the computations and the theoretical predictions. Furthermore, K&L note that a 
pressure gradient of no more than a 0.8 percent per one meter may have been present in 
the experiments. From linear theory, stability results are known to be very sensitive to 
pressure gradients. If a small, unmeasurable adverse pressure gradient was present, then the 
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fundamental wave would reach larger amplitude levels, as recorded. Individually, the effects 
of pressure gradient and frequency differences can not explain the differences between the 
computational and experimental results. So to determine the combined effect of frequency 
and pressure-gradient variations on the fundamental wave F and harmonics 2F, 3 F, PSE 
theory is used to re-compute the subharmonic instability case. 

Using linear stability theory as a guide, a frequency F = 140 and a small adverse 
pressure gradient f3h = —0.02 are selected for another computation using PSE theory. The 
pressure gradient parameter f3h = —0.02 which corresponds to dp/dx = 0.9 percent per unit 
length is of the same order as the possible error given in the experiments. Figures 10 and 
11 compare these new results with the previous DNS results and experiments. As shown in 
figure 10, with this modification of the frequency and base flow, the growth and decay rates 
of the fundamental wave are in much better agreement with the experiments downstream 
of Rs* = 900. The growth rates prior to Rg* = 900 amplify different frequencies in the 
experiments and computations. In figure 11, the computed primary mode (F) and higher 
modes (2 F, 3 F) more closely follow the trend toward decay exhibited in the experiments; 
however, the growth rates of the subharmonic-based modes (F/2 and 3F/2) have increased 
notably as a result of the higher frequency. The earlier PSE and DNS results more closely 
match the subharmonic and its harmonics in the experiments. Although it is unlikely that 
such a large frequency discrepancy from F = 124 noted in the experiments occurred, it is clear 
that the available experimental data is not sufficient to reproduce the results numerically. 

Oblique- Wave Breakdown. 

The final test case used in the evaluation of PSE theory is oblique- wave breakdown. This 
oblique-wave breakdown procedure is due to the nonlinear interactions of a pair of oblique 
waves. Because of this nonlinearity, no adequate formal theory is available to explain the 
breakdown process; however, similar mechanisms have been studied by Hall and Smith (1991) 
with asymptotic methods. Hall and Smith discussed the vortex-wave interactions within a 
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large wave number and Reynolds number limit. To quantify the mechanisms of interest in the 
finite Reynolds number range, DNS and possibly PSE theory are the only options available 
to study the wave interactions. 

Because this alternative route to, or mechanism of, transition is nonlinear, limited re- 
search has been done for oblique-wave breakdown. It is worth mentioning a few of the 
interesting papers that are available. Schmid and Henningson (1992a) studied bypass tran- 
sition by introducing a pair of large-amplitude oblique waves into channel flow. The evolu- 
tion of disturbances was computed with temporal DNS. They found that the development 
of the oblique waves was dominated by a preferred spreading of the energy spectra into 
low streamwise wave numbers, which led to the rapid development of streamwise-elongated 
structures. Schmid and Henningson (1992b) also looked at small-amplitude wave pairs over 
a variety of parameters. They suggest that the mechanism of energy transfer is primarily 
linear. Fasel (1991), Fasel and Thumm (1991), and Bestek et al. (1992) computed this new 
breakdown structure in a compressible boundary layer. They describe the physical structure 
as “honeycomb-like” to identify a distinction from the secondary instability A-like struc- 
tures. Chang and Malik (1992) used PSE theory to examine the breakdown of supersonic 
boundary layers because the dominant first mode is an oblique wave in supersonic flows. 
Chang and Malik found that even waves with amplitudes as small as 0.001 percent that are 
initiated at the lower branch can lead to transition in this breakdown scenario, depending 
on the frequency of the induced oblique waves. Finally, Joslin et al. (1992b) presented 
the preliminary results (computed by both spatial DNS and PSE theory) of oblique-wave 
breakdown in incompressible boundary layers. This paper presents a more comprehensive 
discussion of oblique- wave breakdown in incompressible boundary layers, and PSE theory is 
again evaluated for accuracy in predicting this nonlinear flow development. 

The profiles for the oblique-wave pair are obtained from linear stability theory for the 
Reynolds number R$* = 900, a frequency u = 0.0774, and spanwise wave numbers f3 = ±0.2. 
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Spatial DNS computations are performed on a grid of 901 uniformly spaced streamwise nodes, 
61 wall-normal collocation points, and 10 symmetric-spanwise modes. In the streamwise 
direction, the outflow boundary is 465 8* from the inflow boundary; the far-field, or free- 
stream, boundary is 75 8* from the wall; and the spanwise boundary consists of a length 
equal to one half of the spanwise wavelength, or \ z /2 = 7r / f3. For the time-marching scheme, 
the disturbance period is divided into 320 time steps. For the PSE computational approach, 
100 wall-normal grid points are used; 7 frequency modes and 7 spanwise modes of series (21) 
are used; and the far-field boundary is 58<5* from the wall. 

The input modes are represented by equation (5), which are truncated to four terms. 
The disturbance forcing consists of modes (1,1) and (1,-1), or (cu,/3) and (cu, — /?) , and their 
complex conjugates (-1,1) and (-1,-1). Theoretically, if these modes self-interact initially, 
then only certain higher modes are likely to be excited (supplied energy). These higher 
modes are: (0,0), (0,2), (2,0), and (2,2), etc. 

In the first test case, the oblique waves each have the small amplitude A° 1 = 0.001. In 
figure 12, the primary disturbance (1,1) and the higher modes that were predicted by PSE 
theory are compared to the DNS results. The comparison shows that the modes are in quan- 
titative agreement. Of the modes that were likely to be excited, all received energy initially. 
The streamwise vorticity component (0,2) grows rapidly because of the self-interaction of 
the oblique-wave pair. All other modes grow more slowly downstream than the streamwise 
vortex, and these other modes contain less energy by orders of magnitude. As a result of 
the rapid growth of the vortex mode (0,2), the oblique waves interact with the vortex which 
leads to an amplified harmonic (1,3). This (1,3) mode gains energy sufficient to overtake 
the other initially excited modes, but insufficient to overtake the oblique waves. As shown 
in figure 12, the vortex modes self-interact to supply energy to the (0,4) mode, which has 
roughly the same growth rate as the (0,2) mode. Although the computations were discon- 
tinued, the disturbances will eventually decay and will not lead to transition because the 
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primary oblique waves decay after they pass the upper branch of the neutral curve, and all 
other modes are decaying or becoming neutrally stable. Similar to the previous test cases, 
PSE theory profiles are compared with the spatial DNS for accuracy. In figure 13, stream- 
wise disturbance velocity profiles are shown for the downstream location Rs* = 1178. PSE 
theory predictions agree very well with the DNS results for all modes shown. Note, PSE 
theory captures the near-wall structure evident in the fundamental wave profile in figure 13, 
in agreement with the DNS results. 

For a final computation, the oblique waves are introduced with larger amplitudes A° 1 = 
0.01. The computed primary disturbance (1,1) and higher modes are shown in figure 14. 
Again, the modes predicted by PSE theory are shown to be in agreement with the DNS 
results. Similar to the small amplitude case, the small wave number modes gain initial energy. 
The vortex mode (0,2) is again the dominant, higher order mode. The self-interaction of 
the wave pairs and the interaction with the streamwise vortex lead to a rapid cascade of 
energy to the other modes. Rather than the meager growth and downstream decay, these 
higher modes now grow with growth-rate characteristics similar to the vortex mode. The 
vortex and harmonics rapidly overtake the introduced waves (1,1) and breakdown occurs. 
At breakdown, the spectrum undergoes a filling so both the DNS and PSE computations 
are underresolved near the downstream end of figure 14. Further evidence that the onset of 
transition from laminar to turbulent flow has begun, the skin-friction curve begins to rise. 
Finally, a comparison of local disturance profiles yields the good agreement between PSE 
and DNS comparable to the low amplitude comparison of figure 13. 
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5. Concluding Remarks 


In this paper, PSE theory results were evaluated for accuracy in predicting convective 
disturbance evolution on a flat plate. PSE theory predictions were compared with spatial 
DNS results for 2-D Tollmien-Schlichting wave propagation, subharmonic breakdown, and 
oblique-wave breakdown. 

For 2-D Tollmien-Schlichting wave propagation, the modes predicted by PSE theory 
were in very good quantitative agreement with the DNS results, except a small discrepancy 
in the mean-flow distortion component was discovered and attributed to far-field boundary- 
condition differences. 

For the test case of subharmonic breakdown, the PSE theory results were in very good 
quantitative agreement with the DNS results for all modes, even the mean-flow distortion 
component. Also the present study supports the PSE and DNS comparison made by Herbert 
(1991) for subharmonic breakdown. 

The present computations were in good qualitative agreement with the experiments of 
Kachanov and Levchenko (1984). In light of some modal discrepancies between computations 
and experiments, a possible explanation was given. By introducing a small, adverse pressure 
gradient in the basic flow and allowing for a small effective frequency variation in the distur- 
bance, the computations predicted modal trends which were similar to the experiments for 
the fundamental disturbance. 

For the complicated test case of oblique-wave breakdown, all modes predicted by PSE 
theory were shown to be in good quantitative agreement with the DNS results, even for the 
mean-flow distortion component. Furthermore, these oblique- wave pairs were shown to self- 
interact to excite a streamwise vortex structure, which agrees with the findings of Schmid 
and Henningson (1992a, 1992b). If the initial wave amplitudes are above a threshhold, 
the interaction of these waves and the vortex can lead to a breakdown that bypasses the 
secondary instability stage. Irrespective of the initial amplitudes, the streamwise vortex 
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mode becomes the dominant, higher-order mode. This dominance is significant because the 
presence of small roughness elements may generate oblique wave packets that can interact 
and lead to the increased presence of streamwise vorticity. 

Relevant to the applied engineering community, is the fact that PSE theory has accu- 
rately predicted the disturbance development for 2-D Tollmien-Schlichting wave propogation, 
subharmonic breakdown, and oblique-wave breakdown on a flat plate. A difference in the 
far-field boundary-condition treatment for the PSE method was identified which may cause 
mean-flow distortion variations; however, all other components appear to be affected very 
little as a result of this far-field treatment. In the near future, PSE theory will be a use- 
ful tool for the engineer, and for certain types of flow problems, the use of PSE theory is 
expected to increase. However, because of the assumptions that underlie PSE theory, its 
application to problems other than those for which it has been adequately tested must be 
done with caution. For example, can PSE theory handle the potentially absolute instability 
of a recirculation bubble? Studies of the kind presented in this paper are necessary to provide 
verification of new theories and to study the unknown complex physics of various flows. 
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Figure 1. Amplitude growth with downstream distance, from two-dimensional wave with 
amplitude A° 0 = 0.0025 (r.m.s.). 
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Figure 2. Streamwise (u) mean-flow distortion and fundamental wave profiles with normal 
distance from the wall at the local Reynolds number R$* = 1519. 
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Figure 3. Normal (v) mean-flow distortion and fundamental wave profiles with normal dis- 
tance from the wall at the local Reynolds number R$* = 1519. 
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Figure 5. Maximum amplitude growth with downstream distance for the primary wave 
[F = (2,0)], subharmonic [F / 2 = (1,1)], and various higher harmonics [2 F = (4,0),3-F = 
(6, 0), 3F /2 = (3, 1 )]. 
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Figure 6. Streamwise disturbance profiles at i lie downstream location where locally Rf,* = 
1067 for subharmonic breakdown. [F = (2,0), 2 F = (4,0),-F/2 = (1,1), 37^/2 = (3,1)] 
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Figure 7. Amplitude growth with downstream distance for the primary wave [F = (2,0)] 
subharmonic [F / 2 = (1,1)], and various higher harmonics [2 F = (4,0),3-F = (6,0), 37^/2 = 
(3,1)] at y = 0.26(5. 






Figure 8. Streamwise disturbance profiles at the downstream location where locally R$* = 
1049 for subharmonic breakdown. [u, n = u„.f] 
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Figure 10. Amplitude growth with downstream distance for the previous primary waves and 
a primary wave computed with a frequency and pressure-gradient varitation at y = 0.26(5. 
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Figure 11. Amplitude growth with downstream distance for the primary wave [F = (2,0)], 
subharmonic \F / 2 = (1,1)], and various higher harmonics [2 F = (4,0),3-F = (6,0), 3F/2 = 
(3,1)] at y = 0.26(5. 


42 








Figure 13. Streamwise disturbance velocity profiles at a downstream location where locally 
Rfi* =1178 for the oblique- wave case with A° 1 = 0.001. 
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